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Hele-Shaw flow at vanishing surface tension is ill-defined. In finite time, the flow develops cusp- 
like singularities. We show that the ill-defined problem admits a weak dispersive solution when 
singularities give rise to a graph of shock waves propagating in the viscous fluid. The graph of 
shocks grows and branches. Velocity and pressure jump across the shock. We formulate a few 
simple physical principles which single out the dispersive solution and interpret shocks as lines of 
decompressed fluid. We also formulate the dispersive weak solution in algebro-geometrical terms as 
an evolution of the Krichever-Boutroux complex curve. We study in detail the most generic (2,3) 
cusp singularity, which gives rise to an elementary branching event. This solution is self-similar and 
expressed in terms of elliptic functions. 



I. INTRODUCTION 



The zero surface tension limit of Hele-Shaw flows [T] describes a planar interface between two incompressible 
and immiscible phases propagating with velocity equal to the density of the harmonic measure of the interface. 
It is a famous problem whose importance goes far beyond applications to fluid dynamics. Originally formulated 
by Henry Darcy in 1856 [2] in relation with groundflow of water through porous soil, it describes a large class of 
two-dimensional growth processes driven by a harmonic fleld (a.k.a Laplacian growth) and has led to a number of 
important developments in mathematics and mathematical physics, from complex analysis to random matrix theory. 

A compact formulation of the problem is: let 7(i) be a simple planar curve - boundary of a simply-connected 
domain D. The curve evolves in time t according to 

Darcy Law : Vn{z) ~ H{z), z G 7(t). (1) 

Here w„ is the velocity of the curve - an outward normal vector from D, and H{z) is the normal gradient of a solution 
of the Dirichlet problem in a fluid domain D with a source (a sink) at a distant location (often at infinity in theoretical 
setup; not shown in Figure [l]): 

H{z) = -dnP{z), (2) 
Ap^O on A Pk=0, ---loglzl. (3) 

We use i/(z)|dz| for the harmonic measure on the curve jit). 

Laplacian growth models [3, 4 are characterized by intricate finger-like unstable patterns [5] featuring finite-time 
singularities [SJ |7] in the form of boundary cusp formation (see e.g. Figures |3] |4]). 

As such the problem is ill-defined: the Darcy law stops making sense, since velocity and pressure gradient diverge 
at a critical time in a critical point, where a cusp-like singularity occurs. 

An important feature of Laplacian growth is its integrability. In this case, integrability means that, if the initial 
interface is an algebraic curve of a given order, it will remain so at all times before the critical time. Insights into 
integrabihty appeared in early papers [S] and then were further developed in more recent studies UHl HH [12 [TS] . 

Dynamics near a critical point (i.e., close to the blow-up time) belongs to the class of non-linear problems for which 
any perturbation is singular. In this case, various regularization schemes typically lead to different results. 

This situation is typical in fluid mechanics and, in fact, similar to singularities appearing in compressible Eulcr flows. 
There, perturbing the Euler equations by a diffusive mechanism (i.e., by viscosity), or by a Hamiltonian mechanism 
(i.e., by dispersion) leads to different physics and different flow patterns (as an example compare the dissipative 
Burgers equation and the Hamiltonian KdV equation, which differ by terms with higher-order derivatives). 

A traditional laboratory set-up (called a Hele-Shaw cell after its inventor |T]) consists of two horizontal plates 
separated by a narrow gap, initially flUed with a viscous liquid, where inviscid liquid is pumped in at a constant rate 
at the center of the cell, pushing the viscous liquid, away to a sink. Below we will refer to the viscous liquid as fluid, 
occupying a domain D. An inviscid liquid occupies a complimentary domain D — C \ D. Both liquids are considered 
to be incompressible. The boundary initially moves according to the Darcy law ([l]|3|, where p is the pressure of the 
incompressible viscous fluid. A typical pattern is seen in Figure [l] 
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It seems that the most relevant mechanism of taming singularities for experiments in Hele-Shaw cell (even for 
a developed pattern as in Figure [l] [14^) is surface tension. Then p in (|3]) along the interface is proportional to 
its curvature. This perturbation stabilizes growing patterns (as seen in Figure [T]) but destabilizes their analytical 
and integrable structure. Also, it is not expected that patterns possess universal features if the surface tension is 
significant. 

The surface tension effects can be diminished in non-Newtonian fluids where a viscous fluid is substituted by a 
polymer solution. In this case, singularities are controlled by certain visco-ellastic properties of the fluid [15j. 

An interesting modern setting consists only of the viscous liquid without the inviscid liquid. In experiments described 
in |16j . a thin layer of viscous liquid is positioned on a horizontal wetted substrate, with the free boundary contracting 
by some suction mechanism. [47 The advantage of this setting is that the height of the layer may vary (and indeed 
does vary along the rim |16p. which effectively makes the two-dimensional flow compressible. Below we assume such 
setting that allows the variation of compressibility, and will keep the name Hele-Shaw flow. 

Another interesting setting occurs when air propagates through granular media [17 . There, the boundary is not 
controlled by its curvature (unlike in fluids with surface tension) , but the density of granular media may vary on the 
boundary. 

The traditional approach in studies of Hele-Shaw flow assumes the following order of limits: the density of the fluid 
is set to a constant first, while the surface tension is kept non-zero. Then the limit of surface tension going to zero is 
taken, but proves to be ill-defined; the solution blows up. 

This requires a different model for Laplacian growth, such that the zero surface tension limit ([2| is well-defined 
and does not lead to a singular solution. Such a view for Laplacian growth is suggested by the setting with a thin 
layer on a wetting substrate [16j , or propagating air through granular media [17^ mentioned above, and also by the 
model of Diffusion-Limited Aggregation (DLA) [TB], or iterative conformal maps algorithm [19 . In the latter case, 
the interface evolves by coalescence (aggregation) of Brownian walkers supplied at a constant rate from a point source 
(infinity). The probability of aggregation at a given point is proportional to the harmonic measure of the existent 
aggregate. In the limit when the area of a particle tends to zero, the problem of aggregation becomes identical to 
the ill-defined zero surface tension Hele-Shaw flow. Aggregation continues beyond this limit, producing notoriously 
intriguing fractal patterns. In this case, the zero-size limit for particles is a singular perturbation. Due to a flnite size 
and irregular shape of particles, taking the limit of aggregation model to the continuous fluid mechanics effectively 
leads to a compressible fluid, as in the two settings mentioned previously. 

The geometric regularization suggested by aggregation models is a primary motivation of this paper (and of preced- 
ing papers [9j [TOl HH [HI [13] ) . We will relax conditions that fluids are incompressible and curl- free at a microscale, but 
will require that liquids are curl-free at a larger scale and then study a singular limit when compressibility vanishes, 
setting surface tension to zero in the first place. 
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Based on this idea, we will construct a weak solution of the Hele-Shaw zero-surface tension flow, where: 

- formation of a cusp singularity is followed by expanding shocks - one dimensional lines through which pressure 
and velocity of the fluid are discontinuous. The density of the fluid (which is constant away from the shock) 
becomes singular at the shock, which can be interpreted as a deficit of fiuid; 

- shock lines form an evolving branching pattern. The geometry of shocks is governed by transcendental equations. 



FIG. 2: A growing and branching shock pattern, with one (left) and two (right) generations of branchings. The bold line along 
the negative i-axis represents a narrow viscous finger (fluid). At this scale, the viscous flnger is vanishingly narrow. 



In a subsequent publication, we will describe connections between weak solutions of the Helc-Shaw problem and a 
few related mathematical problems. We will show that: 

- the lines of discontinuities (shocks) correspond to accumulation of zeros of bi-orthogonal polynomials of large 
orders; 

- the shock fronts are anti-Stokes lines for the isomonodromy problem naturally related to the Whitham averaging 
of solutions of Painleve equations corresponding to the integrable models associated with particular types of 
cusps. 

Shock fronts typically form a growing, branching tree, as in Figure [2] In this paper, we provide a detailed analysis 
of the shock pattern solution representing the birth of the first branching event (Figure [2] above) of what will become 
a complicated degree-two tree. This solution is self-similar and represents the local branching structure of a developed 
tree. 

A comment is in order: typically, shock waves in hydrodynamics are associated with phenomena occurring at high 
Reynolds numbers, due to the inertial term in Euler equations |20j . Here we encounter a new situation when a 
discontinuous solution occurs in a viscous flow where the Reynolds number is small and inertia terms are neglected 
in the Navier-Stokes equations (see the next section). For a discontinuity of this kind, perhaps "crack" rather than 
"shock" may be a better name. 

The paper consists of three parts: first we present a few (necessary) known facts about singularities in the Hele- 
Shaw fiow and formulate the problem in terms of inverse balayage - an important concept giving insights to the weak 
solution. Then we formulate the weak dispersive solution which allows shocks, and describe their hydrodynamics (the 
Rankine-Hugoniot condition). Next we formulate the flow in terms of evolution of a complex curve. This formulation 
is the most suitable for computations. Finally, we give a detailed analysis of a generic (2,3) - cusp singularity and 
describe the elementary branching event. This particular solution reflects the nature of a more general weak solution. 

We start by formulating the differential and weak forms of zero surface tension Hele-Shaw flow. 
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II. DIFFERENTIAL AND WEAK FORMS OF HELE-SHAW FLOW 
A. Differential form of Hele-Shaw flow 

A thin layer of viscous fluid of height b occupies an open domain D with a smooth simply-connected boundary ^(t). 
The fluid is drained from a point z = oo at a constant rate Q, Figure [T] 

We assume a quasi-stationary flow and therefore neglect the term ^ in the Navier-Stokes equation 
p (§f + V • Vv — i/V^v) = — Vp, where v is the kinematic viscosity, and p the density of the liquid. We also as- 
sume that the Reynolds number is small Re = |v|5/pi^ — > 0, and further neglect the inertial term there. We obtain 
piyW^v — Vp. Further assuming a PoisseuUe profile for the flow, we replace V^v by its average ~(12/6^)v and neglect 
(V^ + V^)v, to obtain the Darcy law 

pv = ~KVp, (4) 

.2 , 

where K — is called hydraulic conductivity. Evidently, this approximation is consistent if the fluid is irrotational 
V X j = , where j = pv, as it is seen from the last formula. If furthermore the density is assumed to be spatially 
uniform, incompressibility of the fluid V • v = and Q imply that pressure is harmonic at constant density: 

Ap = 0. (5) 

We will relax the condition of constant density when defining weak solutions. The fluid is sucked out at a constant 
rate Q = j x dl, where 7 is oriented counter-clockwise. According to the Darcy law Q, Q equals —K§^ Vp x dl. 
Therefore, at the drain: 

" • ^ log^+O f-) as z^oo . (6) 



2tiK r{t) \z 

If the boundary is smooth, the pressure on the boundary is controlled by the surface tension p = a x curvature. If 
surface tension vanishes, the pressure is constant along the boundary p = 0. Then the pressure is a solution of the 
Dirichlet problem in the fluid. 

The time dependent function r{t) (known as conformal radius), or capacity C{t) = r{t) [21 , where i? is a cell 
radius, are important characteristics of the flow. They monotonically grow and have the following hydrodynamic 
interpretation |22j : let us compute the power required to drain the fluid: Nit) = — Vp • j dxdy. Darcy law yields 
that the power is the Dirichlet integral N{t) = K f^{Vp)^dxdy evaluated as N{t) = Q^/{2'KK){logR - logC{t)). 
Then the power is a monotonically decreasing function of time. 

We will see that, as a flow approaches a (2,3)-cusp singularity and enters into a shock regime, capacity increases 
non- analytically, as a square root of time measured from the critical point. After a singularity, capacity again increases 
as a square root of time measured from the critical point. A notable fact is that the coefficient of the square root 



singularity changes discontinuously before and after the transition (see Sec. VI D I. Physically, it means that the power 
required to produce the flow goes through a discontinuous change every time the shock front branches, such that more 
power is required than it would be in the absence of singularity formation j22j . 

B. Weak form of Hele-Shaw flow 

Once the flow reaches a cusp singularity, it cannot be continued any longer. The Darcy law written in differential 
form loses it meaning at the cusp - harmonic measure there diverges. 

This situation is typical for differential equations |23j with conservation laws of hyperbolic type: 

dtu + dJiu) = 0. (7) 

Cauchy problems with such conservation laws are ill-posed, and smooth initial data evolve into singularities devel- 
oping in finite time (such as shocks, vortices, etc.). Darcy equation is of this type. 

The origin of this phenomenon is that conservation equations of hyperbolic type are approximations of well-defined 
problems. A deformation of these equations by terms with higher gradients, controlled by a small parameter h, prevents 
formation of singularities. If a smooth solution of a deformed equation un leads to a space-time discontinuous function 
u\ti^Q as the deformation is removed, it is called a weak solution. A discontinuity (a shock front) travels with a velocity 
given by the Rankine-Hugoniot condition - a weak form of the differential equation ([t]) [201 [H] : 

V=pl. (8) 
disc u 
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In most cases, simple physical principles determine a weak solution without actually specifying a regularizing 
deformation. The best known example is the Maxwell rule determining the position of shock front |20j . or more 
generally, the Lax-Oleynik entropy condition for a viscous solution of equations with hyperbolic conservation laws 

m- 

Here we assume the same strategy. We will be looking for a weak solution of the Hele-Shaw problem when the 
Darcy law is applied everywhere in the fluid except a moving, growing and branching graph T(t) <Z D of shocks (or 
cracks) , and where pressure suffers a finite discontinuity - a weak form of singularity. 

Below we formulate a few simple, natural physical principles which will guide us to obtaining a unique weak solution, 
that we call a dispersive solution. 

Shocks may have different physical meanings depending on the various experimental settings discussed in the 
Introduction. In the two- liquids setting (traditional Hele-Shaw cell), shocks are narrow, extended channels where the 
inviscid fluid is compressed and carries vorticity. In the single fluid setting, shocks are areas with a deficit of fiuid. 
Shocks communicate with the bulk of the fluid by supplying/removing fluid to/from the bulk. Obviously, for shocks 
to occur, the uniform fluid density condition must be relaxed. A difference between shocks and a regular boundary 
is that both pressure and velocity have steep gradients across a shock. Also, the pressure is generally not a constant 
along a shock line. 

In the following, we summarize the weak form of Hele-Shaw problem, postponing the details to Sec. |V] The weak 
form reads: 



the differential Darcy law (^|6l holds everywhere, except on the graph of shocks, where pressure jumps and the 



density (a constant everywhere else) has a single layer density deficit p(x,t) = po — 6{x;T)a{x,t) with a hne density: 
cr(x, t) > (here (5(x; F) is the delta-function on a shock). Both pressure and line density vanish at shock's endpoints 
efc as -y/jx — Gfej. Shocks move with velocity whose normal component Vj^ is directed towards higher pressure regions, 
and obeys the Rankine-Hugoniot condition: 

cr|V_L| = if|discp|. (9) 

We discuss the motivation, interpretation and consequences of this solution in sections below. The conditions imposed 
are restrictive. They determine a shock pattern with a given number of legs, up to a finite number of deformation 
parameters. Shock graphs with one and two branching generations are represented in Figure |2] 



III. FINITE-TIME SINGULARITIES IN HELE-SHAW FLOW 

The dynamics described by the Darcy law ([l] [2| is ill-defined. This has been understood from different angles 
m |5S], and can be seen from a simple argument: let w{z) be a conformal univalent map of the fluid domain D 
onto the exterior of the unit circle, \w\ > 1; the solution of the exterior Dirichlet problem then reads 

p{z) = -hog \w{z)\, v=^\w'{z)\, zel). (10) 

(From now on we set K — l,Q — tt. In these units the area is equal to tt x time. In this section we also set p — 1). In 
other words, velocity of a boundary point is proportional to the density of harmonic measure i/(z)|dz| — Iw' {z)\\dz\ 
of that point. As a consequence, a sharper (more curved) part of the boundary moves with a higher velocity than the 
rest, such that the sharper part is getting even sharper, until it becomes singular. This occurs for generic initial data, 
although a number of exceptions are known [26^ . We will demonstrate this mechanism on a simple but representative 
example, but first we remind the major feature of the dynamics [5][5S]: let us denote by 

t, = -lJ^z-''d.dy=^J^z-Hdz, (11) 

the harmonic moments of the fluid, and let the area of the domain D not occupied by fluid be irto / p = Q-t. Parameters 
{tk} are called deformation parameters. It follows from the Darcy law that [S] 

ik = 0. (12) 

Let us write the inverse (time-dependent) map as 

z{w) = rw + ^UkW^'', |w| > 1, (13) 

k>l 
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where the coefRcient r, chosen to be real positive is the conformal radius ([6|. According to a standard formula of 
complex analysis, the coefficients of the inverse map determine the area of the domain D, which we denote by nt, 
through the area formula [2T] 



(14) 



fc>0 



Let us consider a simple example of a domain D whose shape is given by a symmetric hypotrochoid - an algebraic 
curve with a three-fold symmetry such that all harmonic moments vanish = 0, except for the third one (see 
Figure [3|: 



z{w) = rw + 



„2 ■ 



(15) 



The ratio between the square of conformal radius r and u can be expressed through t^, a,s u — St^r^. Then the area 
formula reads t = — 2(3|i3|)^r'*. 

Clearly, this polynomial in reaches the maximum tc — r'^/2 at rc — l/(6|i3|), and cannot increase further. At 
this moment, the hypotrochoid develops three simultaneous cusp-like singularities. Figure [3j 

The mechanism of the evolution toward a critical point is: a critical point of the inverse map z'(wc) = is located 
outside of the fluid at Wc = {2u/ry^^ = {r{t) /rcY^'^ < 1 and moves in time towards the boundary. At a critical time, 
the critical point reaches the boundary, at w = 1. 

In the critical regime, when the critical point of the conformal map is already close to the boundary z(w)\^^f,'<t' ~ 



2^ — x{(f>) + iy{(j)), the boundary is locally approximated by a real degenerate elliptic curve: 



X(0) 



2 X 

3 rr 



= -4{X-e{t)) [X + 



e{t)^ 



y 



2f 



^(1 



r) 3 



3e(t)0, 
2 



1 - 



1/2 



2 ) ' 3 
At the critical point t = t^,, the curve degenerates further to a (2,3) cusp: 



Close to the critical point, the conformal radius depends on time in a singular manner as r/r^ — 1 
so that the time dependence of power becomes non-analytic: N — ^ {1 — t/tcY^^ . 



(16) 
(17) 

(18) 

-^-{l-t/tcf\ 



Once the (2,3) singularity occurs, the shape of the pre-cusp finger is given by the elliptic curve ( 17l. Up to a scaling 



of coordinates and time the curve is universal - it does not depend on details of the original domain: 



YiX,t) 



1 

tr. 



3/4 



Y 



X 



(19) 



Self-similarity is preserved after the flow passes through the critical point, t > tc- 

This simple analysis has a straightforward generalization for all cases where a finite subset of {tk}k>o is non- 
vanishing. At least for this class of initial conditions, the flow leads to cusp-like finite-time singularities. For purposes 
of analysis of singularities it is sufficient to consider only this set of domains. Of course, the class of initial conditions 
which leads to singularities is much wider. We do not attempt to classify it here. We only mention that domains 
with a flnitc number of non-zero exterior harmonic moments belong to the class of generalized quadrature domains. 
For generic properties of generalized quadrature domains, see e.g., |27j . 



A. The first classification of singularities 



The fact that for rather general initial data, the boundary evolves to a cusp-like singularity had been noticed already 
in the earlier papers [25|, I28j where the model of Laplacian Growth was originally formulated. Further developments 
are found in [51 [3], and [31 El UHl EHl Ell El] • For the most generic initial conditions, a cusp (2,3) occurs. Higher-order 
cusps and corner-like singularities require special initial conditions. 

In |11| . it was argued that, if a complement of a fluid domain £) is a generalized quadrature domain, then cusps of 
the type (2, 21 + 1) : ~ are possible, and that in the critical regime a finger developing into these cusps is 



FIG. 3: Symmetric hypotrochoid ( 15 1 evolving under Darcy law (left) reaches the (2,3)-cusp singularities (right), when all three 
critical points (red dots) hit the boundary at the same time. The shaded contour lines are the equi-pressure lines. The dashed 
lines are the stream lines. The red dashed lines are branch cuts -a skeleton. 




FIG. 4: Deformed hypertrochoid z{w) with Ar — r — Ar'^ — y ^ — 8r* + IGr'', where one critical point 

gives rise to a (2,3) singularity while two other critical points remain in the domain D. Closed lines inside the critical boundary 
show the pre-critical evolution. 



described by a real degenerate hyperelliptic curve 

I 

Y'^-4{X^e{t))l[{X~d,{t))'. (20) 

1=1 

This means that a real oval of the curve - a set where both X and Y are real - is a graph of the finger {x, y) — {X, Y). 
The double roots di (double points) are all located in the fluid domain and simple critical points e and — oo are located 
outside of the fluid. 

The condition that no critical points are found in the fluid, necessarily means that the exterior critical points 
coincide, giving double points, such that the complex curve is degenerate. A cusp of type (2, 21 + 1) occurs when the 
branch point (real root e) merges with I double points. 



For reference, we give explicit formulas for the generic hyperelliptic curve (20 1 following the Hele-Shaw evolution 
[llj . The curve depends on / deformation parameters {t2k+i},k = I,..., I (not to be confused with the harmonic 
moments tk), and it is uniformized by the formulas 

X = e(^)-0^ Y^E':=\in+h)t2n+iEVo (fe^ ^ ^ (21) 
h + E[t\ ^^i2.+i(-e(t))'= = 0. (22) 
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Here ti is a time before the cusp singularity, ti ^ t — te- 
la. [29] it has been argued that dynamics can be continued through a cusp of the type (2, 4A: + 1) for the price of an 

appearance of an unstable point when a boundary touches itself. Another stable solution was found in [T3] when the 

fluid undergoes a genus transition: it becomes multiply connected. 

However, in the case of the most generic cusp (2,3) and its descendants (2,4/c — 1), fc > 0, no smooth solutions are 

possible. This situation is the subject of this paper. 



IV. ALGEBRO-GEOMETRIC FORMULATION OF HELE-SHAW PROBLEM 



A. Real algebraic curve 

We briefly review some elements of the algebro-geometric description of the Hele-Shaw flow. It is based on the 
notion of real complex curve. 

A boundary of an algebraic domain (the only case we consider) is described by an equation /(z, z) = 0, where 
/(z, z) = X)nm o.mn.z^z™' IS a polynomial. It is proved to be convenient to consider a Riemann surface or real 
algebraic curve defined by the equation /(z, z) = f(z, z) — 0, where z and z are two complex variables. For example, 
a curve for a symmetric hypotrochoid (3) is /(z, z) — {zzY — 4:rc{z^ + z'^) + 4(r^ — ^)(1 + ^){zz) — 4r^(r^ — ^) 

Among many solutions of the equation /(z,z) — with respect to z (sheets of a curve), there is one z = S{z) (a 
physical sheet) which describes a boundary where z = z, 

S{z) = z, z e 7. (23) 

It follows from Darcy law that the time derivative of function S{z) is holomorphic in the fluid. On the boundary 
and in the fluid domain D a complex velocity v — Vx — iVy is obtained by 2v — dtz — dtS{z). Everywhere in the fluid 
the velocity is holomorphic. In terms of the function S the Darcy law reads: 

2v{z,t)=dtS{z,t), (24) 

dtS{z,t)^2id,(^{z,t), (25) 

where 

<P{z,t)^il;{z,t) + ip{z,t) (26) 

is the potential of the flow, whose imaginary part is the usual scalar pressure p, and the real part is the stream function 
"0. Wherever pressure is harmonic, the potential is an analytic function. 

As we have said, to study local behavior of an isolated cusp, it is sufficient to consider domains with a finite set of 



non-zero harmonic moments (11 1. In this case the function S has a multiple pole at infinity. It is represented by a 
truncated Laurent series with respect to a drain at infinity and an arbitrary regular point, say z = 0, outside of the 
fiuid: 

00 

S{z)^ S+{z) + S-{z), 5_(z) = ^z-l + ^^;feZ-^ (27) 

k>l 

where we set the positive part to be a polynomial: 

5+(z) = ^fctfcz'=-i. (28) 

fe>0 



It encodes information about the moments of the exterior [3T and from (24 1 it docs not evolve in time (12 1 



S+{z,t) = Q. (29) 

Parameters of the flow - time t and deformation parameters tk - can be seen as residues of the differential z~'^S'(z)dz 
at z = cxD (a drain) . 

Before a critical time, complex velocity of the fluid is holomorphic everywhere in the fluid domain. In this case, 
the negative part of the function S{z) is an analytic continuation of z to the fluid domain: 



27ri 



D 
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In this case, S{z) is called a Schwarz function f3P. All singularities of S-{z) are branch points located outside of the 
fluid domain. 

A generating function ~ integral of the meromorphic differential df2 = S(z)dz from a point of the boundary of the 
fluid Zry to a point of the fluid 



n{z) 



S{z')dz' 



(31) 



is a convenient way to describe Hele-Shaw flows. For example, the boundary of the fluid given by (23) is a solution 
of the equation 



|z|2 2Ren{z). 

A graph of —\z\'^ + 2Keil{z) close to a critical point is plotted in Figure [s] 



(32) 



FIG. 5: The plot of —\z\'^ + 2Re / S{z)dz before a critical time (left) and at a critical time (right). The height difference has 
been blown-up by scaling the vertical axis and by applying arctan. 



In a critical regime, i.e., close to a cusp, it is convenient to describe the boundary in Cartesian coordinates Y = 
5j(z — S{z)) and X = 5(2 + S{z)), and treat F as a function of X. In these coordinates pTj Darcy law (24l reads: 



dtY{X,t) ^ -dx<l>iX,t). 



(33) 



This follows from an important property of the Darcy law (24 1: the form doj = Sdz + 2i(f>dt is closed. Therefore, the 
Darcy law is invariant under canonical transformations {z,S) — > {X, Y) : duj = —2iYdX + 2i(f)dt, up to an exact 
form. 

In a critical regime, it is also convenient to choose a critical point as origin and to redefine time t ^ ti ^ t — tc 
and deformation parameters as t2k+i = 2(2^+1) ^'*^^°° {-X)-^-'^/^YdX. Then, in the case of a (2,2^ + 1) - cusp, a 
hyperelliptic curve (20 1 is a Laurent series with respect to the local parameter at infinity X^/"^: 

Y{X) = Y+{X)+Y_{X), 
Y+{X) = T!u2,{k+\)t2U+i{-Xf-y\ 
^(X) = -(-X)i/2 + g(-X)-i/2 + 0_(X), 

where Y- and 0_ consist of further negative powers in X^/^. Here, C is the capacity of a finger ~ relative to the 
capacity at a critical time (capacity for a compact set is defined just below ([6])). Darcy law then states that the 
positive part is conserved: 



(34) 
(35) 
(36) 



y+(X)=0. 



(37) 



Eqs. (20 1 give solutions of the Darcy law written in the form (33 1 under assumption ( |35[ ). The parameter 
(treated as a function of X) is the complex potential of the fiow at the point z = X + iY{X). 



in (211 
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B. Skeleton and inverse balayage 



If the number of non-zero harmonic moments is finite, the function S{z) (27 1 has only a finite-order multiple pole 
at infinity and branch points in the finite part of the complex plane. Before a critical point forms, all singularities are 
moving branch points located outside of the fluid. Their dynamics encode the entire flow. 

In this case, there is a unique, special way to draw branch cuts. They can be drawn as a (possibly multi-component) 
curved graph such that the jump of the differential — idf2 = — iS'dz, being canonically (counterclockwise) oriented 
along the curve through every branch cut, will be real and positive: 

-idisc S'(z,t)dz = 2crfc(z,t)|dz| > 0, zeY.k{t), (38) 

where k labels branches of the graph. Here the differential is taken along the cut and \dz\ is an element of an arc-length 
along a cut. The canonical orientation (important for all signs in formulas below) is depicted in Figure [6] 



If S is treated as a complex vector, then (38) reads 



disc S ^ 2CTfen. (39) 



disc n = Jp S{z) dz 

disc S = ^(iip) — 5(down) = 2\a 



ij ■ Idisc 51 , • 

I = 1^. — = m 



X- 
down 



FIG. 6: In the left panel, we show the three possible stages of evolution: i) when a branch point is outside the fluid, ii) when a 
finger forms a cusp (a branch point is on the boundary) , and iii) when a shock develops after a cusp (branch points are inside 
the fluid). Canonical orientation of the contour (dashed line F) define how the discontinuities are taken. In the right panel, 
we show the normal velocity of the shock by the thick arrows and the flow of fluid around the shock by dashed arrows. The 
arrows illustrate that vorticity carried by a shock is compensated by the circulation of the fluid around the shock. 



This fact is equivalent with the following analytical continuation of S'_ (z) (|30p to the domain D 



S-{z) = 



D 



E 



(40) 



with |dC| the arclength measure. In other words, 5_(z) is the electric field produced by a single layer of charges 
CTfc distributed along a graph E. An important theorem on generalized quadrature domains |36] tells that u exists 
and is unique. The curved and generally multi-component graph S, confined in the domain D, is called skeleton or 
mother-body [36 . Line densities of the skeleton depend on time only through time dependence of the branch points. If 
a branch point e{t) is simple (which we assume), the line density a{z,t) vanishes at the branch point by a semicircle 
law as 



a{z,t) ^ V\^~e{t)\ 



(41) 



In the case of the hypotrochoid (15), the graph consists of three curved lines connecting the branch points. They 
are drawn as in Figure |3] In the case of symmetric hypotrochoid (15), the branch cuts are three straight seg- 
ments, connecting the branch points e(t) = {r(t)/rc)^^^ x (1, e**^'^/"^), with line densities at the branch points 
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4(2/3)-^/^rc(l — r^/r^)(rc/r)^/'^[^|2; — e(t)\ + (4/3)(rc/r)''/^|z|/rc-\/|2; — e(t)|] (where z is measured from a position of 
a cusp). 

The major property of generalized quadrature domains (actually serving as definition) |27| is: area averaging of 
an integrable analytic function /(z) over any domain B is reduced to counting of singularities of the function S{z) 



h Ib /(-^)^^'^ — 2^ ias This follows from the representation (40 1 



^/B/(^)pdxdy = :^Jggf{z)S{z)dz 

(42) 

= /{S}ni3^(^)/(^)|d^l- 

If the domain B does not contain moving singularities or the drain, then the average stays constant in Hele-Shaw 
flow. Conversely, if the domain B contains moving singularities, or a drain, the average depends on time. 

The procedure of sweeping a domain B into a line graph {E}n-B with a line density a is known as (inverse) halayage 
[55] . A Newton potential created by uniformly distributed mass outside of the fluid, measured inside the fluid, in D, 
will be the same as the Newton potential created by the (non-uniform) line densities a. 

The balayage procedure allows the following insightful interpretation of Hele-Shaw flow. The "skeleton" - moving 
branch cuts - may be considered as time-dependent sources of the viscous fluid. In this interpretation, the fluid has 
no boundary, only time-dependent line sources. An extended fluid occupies all the plane except the sources and moves 
with velocity v — {\/2)dtS. The pressure of the extended fluid is deflned hy v = — Vp . The original boundary can 
then be restored (if necessary) as a real oval of the complex curve - a planar curve where S{z) = z, or (equivalently) 



a zero-level set of the pressure p = 0. Darcy Law (24 25 1 connecting velocity and a complex potential then applies 
both inside and outside the fluid. 

In this interpretation, the entire flow is encoded by the motion of the skeleton, or rather by the motion of end 
points of the skeleton. Velocity, pressure and stream function generally have finite discontinuities across the skeleton. 
Also, pressure may not be constant along the sides of the skeleton, while velocity is not normal to the skeleton. A 
singularity occurs when a growing skeleton intersects the boundary of the fluid. 

In the next section, we formulate the weak solution of the Hele-Shaw flow and derive equations for the moving shock 
graph. The same equations are applied for a moving skeleton. Vice-versa, an evolution through a cusp singularity 
will be "continuous" only if shocks are governed by the same law as the skeleton. 

V. WEAK SOLUTIONS OF HELE-SHAW PROBLEM 

A. Weak form of Darcy law: the Rankine-Hugoniot and admissibility conditions 

Now we assume that after a singularity has occurred, pressure is a harmonic function everywhere except on a moving 
shock graph T(t) C D located in the fluid domain, where pressure and stream function have finite discontinuities. 



Then so does the r.h.s. in (25 1. 

The function S{z) = S-\-{z) + S-{z) that we defined before the cusp occurs changes its property and is no longer 
a Schwarz function. In particular, it exhibits a discontinuity in S-{z) along the shock (in addition to the skeleton) 
such that: 

-idiscS{z,t)dz = 2ak{z,t)\dz\>0, zeTkit). (43) 



We keep the same notations for line densities along shocks as lines densities along a skeleton ( 38 ) 



We also assume that shocks that appear in the post-critical regime do not affect the fluid far away from the critical 



point. We require that the lines densities are real and positive, as for a skeleton (Sec. IV B I. This is the major 
condition. We illustrate its meaning below. 

Therefore, after the critical time, S{z) develops branch cuts in the fluid domain D in addition to the branch cuts 
in D. We write: 

With an assumed orientation as in Figure[6] we can interpret (44) as a deficit of the fluid density: p = pa — S{z;T)a{z) 



A comment is in order. Before a critical time, the function S{z,t) in (|43| was analytic in a vicinity of the fluid 
boundary, and therefore was a Schwarz function of that boundary. Not so anymore after a critical point. It jumps at 
branch cuts that protrude from a skeleton. For the same reason, a domain D (complement to the fluid domain which 
includes shocks) is not a a generalized quadrature domain any longer. 
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Let us integrate (25) over a closed loop dB bounding a domain B drawn anywhere in the fluid. If the loop does 



not cross a shock, we obtain 



Ini 

2dt 



1 d_ 
2di 



Re 



S{z)dz 
S{z)dz 



I (V-j) d^z= i j xdl= </dV, 

JB JdB J 



(V X j) d^z - 



jdl 



SB 



dp, 



(45) 
(46) 



where V x j = dyjx — dxjy is the vorticity field, and V • j = dxjx + dyjy is the divergence of velocity field {2dj — 
V • j + iV X j), and j = pv. The first two equalities in each line of (45 46 1 are identities, the last is one is the Darcy 
equation. 

Let us now consider the integral S{z)dz over a loop which does intersect a shock graph at one or two points. 
This integral measures the density of the portion of a shock surrounded by the loop and under (43 1 stays purely 
imaginary at all times. We conclude that 



dt 



Re (p Siz)dz = 0. 



(47) 



Letting the contour now shrink to an infinitesimal loop, we obtain the differential form of (47 1 as: 



— Re [discs' dz]r = 0. 



(48) 



The total time derivative in (48 1 has two contributions: one from the time evolution of the function S{z), the other 



from the motion of shocks. We denote the velocity of the shock front (normal to the instantaneous curves F), by Vj_, 
directed along the vector n as in Figure |6] Then the total time derivative ( [is] ) becomes 



Re [discs dz + V II (discs' •V_L)|d.z|] r = 0, 



(49) 



where Vy represents the derivative along the direction tangent to the front along the vector {. as in Figure |6] 

For the first term, we use a kinematic identity Re [disc S dz] = 2disc wy |dz|, valid on both sides of the shock. From 



the reality of the jump (43 1 it follows that the second term in (49) is purely real, and equals 2Vy (crVj,) |dz|. 
Together, it yields to flie condition 

V|| Jl + disc jii =0, J_L = crVj_, jii = pqvh. 



(50) 



The first term in this equation represents the transport of mass due to motion of the shock (normal to the shock 
itself), while the second is the vorticity of the surrounding fiuid fiow. They compensate each other. This condition, 
derived solely from the requirement that a is real, suggests to interpret a shock as a single layer of positive vortices. 
Then a is the (smooth part of) density of vortices. The core of the vortex is of the order of h^^^ ~ the width of the 
shock. Then the reality condition for a means a zero -vorticity condition for the fiuid. 

Using Darcy law we replace the fiuid velocity in (50) by — Vyp, and integrate (50) along the cut. We obtain the 
Rankine-Hugoniot condition 



crVj^ = (disc p) n. 



(51) 



(the constant of integration is fixed from the assumption that both the line density and discontinuity of pressure 
vanish at endpoints of the cut). The fact that cr > is positive means that a shock moves toward the direction of 
larger pressure. 



Furthermore, calculating the discontinuities on both sides of the Darcy law and using (43) we obtain 



(7 = —disc V w'tp 



(52) 



This formula has the following interpretation. Let us assume the shock filled with some material (say, an inviscid 
fiuid in the Hele-Shaw cell), of line density a. Then the continuity condition (T + V||(aV||) = 0, tells that this material 
is moving along the shock with velocity V||. Combining this with (52) and integrating along the shock, we obtain a 
counterpart of ( |50[ ) : divergency of the sliding current along a shock equals the discontinuity of current of fiuid normal 
to the shock (see Figure [6]): 



V|| J|| -disc j_L = 0, J||=o-V|| 



(53) 
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Integrating this formula along a shock gives a counterpart of the Rankine-Hugoniot condition: 

ctV|| = (discV)^- (54) 
At an end point of the shock e{t), the total velocity of the shock matches the velocity of the branch point: 

e(i)- (V^+V||)U^,(,). (55) 
We notice that (at least close to the endpoint), the stream of mass is directed towards the branch point. 



The Rankine-Hugoniot conditions ( |5l[ |54[ ), inequality (43 1 and the differential form of the Darcy equation (251 
combined give the weak forms of the Darcy law. 

The same conditions hold for a skeleton before and after critical time. After a critical time, shocks can be considered 
as continuation of growth of a skeleton inside the fluid. From this point of view, skeleton and shocks together form a 
graph which grows under the same law. As we will see below, the boundary is singled out as a position of the flrst 
branching event. Following branching events already occur in the fluid as branching shocks. 

Summing up: three conditions 

a'V± = (disc p) n, a = —disc V|| V', cr > 0, (56) 

together with the differential Darcy law outside of the shock's graph constitute the weak form of Darcy law. 

We note that this weak solution only applies to algebraic solutions unlike the conventional Hele-Shaw problem. 



Though condition (56 1 tells one how the shocks and domains evolve in the next moment, the determination of the 
driving force — pressure — is not given by the simple Dirichlet problem as in the conventional Hele-Shaw problem. To 
obtain pressure one needs to know S{z), which contains a finite number of branch cuts and poles. If the function S{z) 
is algebraic at some moment in time, integrability of the problem guarantees that it will remain so at all times. 



B. Algebro-geometrical formulation of the weak solution 

Having real complex curve in mind (not only its real oval) has proved to be useful for studying a Hele-Shaw flow 
and is necessary to understand its weak solutions. A Hele-Shaw flow of an algebraic domain is then understood as 
the evolution of a meromorphic differential dil = S{z)dz. The properties described above can be cast as algebro- 
geometrical properties of an algebraic curve: 

1. the complex curve is real; 

2. the condition that all densities are real yields 

Re f dfi = 0, for all cycles in the fluid; (57) 



3. shocks form a graph determined by the condition 

Rediscr2 = 0, (58) 

and the admissibility condition |4] 

4. the condition that all densities are positive yields that Re f2 is increasing away from the cut; 

5. the meromorphic differential dil has multiple poles at infinity, the residues at infinity of the differential z~'^d^l 



are the deformation parameters (11), and time is the residue of the differential dfi at infinity. 



The graph (58 1 consists of curved lines, sometimes called anti-Stokes lines. Condition 3 selects admissible anti-Stokes 
lines. 

Complex curves obeying condition 2 are called Krichever-Boutroux curves. Evolution of a curve obeying the 
Krichever-Boutroux condition through a change of genus had already appeared in studies of Whitham averaging in 
integrable systems by Krichever |37j . Similar condition has been introduced by David in the theory of 2D quantum 
gravity |38j , and recently was recognized in studies of distribution of zeros in orthogonal polynomials by Bertola and 
Mo 131]. Neither appearance is a coincidence. Semiclassical asymptotes of orthogonal polynomials and Whitham 
genus changing transition are ultimately related to the Hele-Shaw flow. 

According to property 3, the real part of the generating function ri(z) = 5'(C)dC, where z^ is a chosen point on 

the boundary, is a single- valued function. The integral of the differential over a cycle involving a drain is f dfl = 2it. 
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Being cast in this form the properties of the curve fully and uniquely determine a smooth Hele-Shaw solution before 
a singularity is reached, and a weak dispersive solution afterwards. 

As we have discussed, the boundary of the fluid (before it reaches a singularity) is a degenerate curve. A singularity 
occurs when a branch point meets a double point on a boundary, forming a triple (or higher order) point, so the curve 
further degenerates, while the boundary of the real oval forms a cusp. Beyond this time, the triply-degenerate point 
splits into three regular branch points, so that the curve becomes non-degenerate. In other words, going through a 
singularity, the curve changes genus. This process is illustrated in Figure [7j 



In a critical regime (i.e., close to a cusp-singularity), a complex curve locally is hyperelliptic = R2i+i{X) (21 1 
and 

fl{X) = -i / YdX, (59) 

^ — oo 

where —oo is a distant point on a skeleton (negative real axis). In this case Im on both sides of the cut has opposite 
signs, while Re is zero on the cut. For hyperelliptic curves, conditions 3 and 4 arc formulated in terms of level lines 
of the field of the potential Re fi: 

Ren{X)>0, X^T; (60) 
Rcr2(X)|r = 0. (61) 

Shocks and skeleton are zero level lines of the potential Re fi, such that the potential stays positive in the neighborhood. 
Then the line density of shocks and the mass deficit accumulated by a shock between two points Xi and X2 are: 

a{X) = \Y{X)\r, (62) 

^2 X, 



a\dX\ = \hnn{X)\ . (63) 



In the next section, we show how these two conditions uniquely determine the evolution of the elliptic curve through 
the critical time. 



VI. BEYOND THE (2,3)-CUSP SINGULARITY: SELF-SIMILAR WEAK SOLUTION 

Here we give a detailed analysis of the generic (2,3)- cusp singularity and illustrate the nature of weak solution. The 
computations below describe the evolution of a unique elliptic Krichever-Boutroux curve. We emphasize that part 
of this analysis appeared in the papers [3H1 SOI HI] , which is devoted to Whitham averaging of the triply-truncated 
solution of Painleve I equation. In this section we again set po = 1. 

We study the evolution of the elliptic curve ( 16p9 1 



y2^_4(X-e(i))(^X+^j , (64) 

through the critical point where e{t)\t=o = 0. It is in order to remind how complex coordinates X and Y are related 
to a coordinate z = x + iy of the fluid. It is z/vc = X + iY{X). On the boundary of the fluid X and Y are real 
and they coincide with Cartesian coordinates {x,y) of the fluid. Elsewhere X is complex. The coordinate of the fluid 
is given by the map X z. The map covers the physical plane twice, as shown in Figure [8] so one must choose a 
physical branch of the map. Under this map a point in the fluid z is in one-to-one correspondence with a point on the 



X-plane, except on the fluid boundary. There, two points z and z correspond to the same X. As it follows from (16 1 
X and Y{X) have different scaling X ~ x/rc, Y ~ (x/rc)'^/^. Therefore at x <^ Vc i.e., at the critical regime where 
approximation of the finger by the elliptic curve make sense, Y{X) <C X and therefore z « X. Having this fact in 
mind below we express the flow in the X-plane. There the boundary of the fluid - a narrow flnger - is represented by 
its skeleton - a cut on the negative axis ImX = 0, Re A" < e = 63, as in Figure [7] 



A. Elliptic curve — genus transition 



Figure [7] illustrates the genus transition of the curve and hydrodynamics of shocks. Before the critical time, two 
branch points coincide in the fluid. There, pressure and velocity are smooth. The branch point 63 is the tip of the 
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FIG. 7: The equi-pressure lines and the flow (dashed) hnes before the critical time (left), and after the critical time (right). 
Pressure of fluid gets larger as the shade gets darker. The boundary of the fluid is a very narrow finger ~ around the 
skeleton line (the thick line on the real axis) and it is not depicted. The red dots are the branch points, 61,2,3- Real 63 is the 
tip of the finger. Before the critical time 62 = ei is the double point (larger red dot in left figure). The thick lines (right panel) 
connecting the branch points are the shocks. They emanate with the angle 27r/3 and diverge at their ends by a slightly large 
angle 0.09184987r + In the figure on the right, the orange and green dashed lines that emanate from the branch points are 
not-admissible level lines of Ref2. The fluid flows to the lighter region, toward low pressure. Shocks move toward the darker 
region (higher pressure). 



finger. At the critical time, the tip meets the double point. After the critical time, the double point splits into 
two branch points. They are the endpoints of the shocks. The tip of the finger 63 retreats. Fluid goes towards 
lower pressure (lighter regions), creating a deficit on shocks where fluid decompresses. Shocks move in the opposite 
direction, towards higher pressure (darker regions). 
We write the curve as: 

y2 ^ _(4^3 _ _ ^ _4(^ _ e3(i))(^ - e2{t)){X - e^{t)). (65) 

where e^{t), e2(i), ei(t) are time-dependent branch points. One root, say 63, can always be chosen to be real. The 
other two are complex conjugated 62 = ei. The coefficient in front of does not depend on time and can be removed 
by translation, such that 63 -|- 62 + ei = 0. 

The positive part of the curve is ±y+ = 2iX^/^ — i^X^'^/'^. The coefficient of the X^^/^-term is is proportional 
to time. We normalize it as — j(?2 = 3t (this choice corresponds to ti = 6t in (351, which further requires a change 
(j) (p/6 below). Under this setting, 

-|g2 = -e2-K|e2p = 3t, .93 = 4e3|e2p, e3 = -2Re(e2). (66) 

1. Before the critical time (t < 0) the real branch point 63 is negative. It is located on the boundary of the fluid. 
Therefore, Re e 1,2 > 0. They lay in the fluid. Condition that the curve must have no branch points in the fluid 
(following from incompressibility of the fluid) requires that two branch points coincide to a real double point: 



62 = ei = —63/2 > 0, so that the curve is degenerate, as in (17 1. Condition (66 1 yields: 



63 = -2V-t , 62 = ei = V-t. (67) 

The double point 62 = ei located in the fluid and the branch point 63 located outside of the fluid are depicted 
in Figure [t] The curve reads: = _(4x3 + l2tX + 8{-tf/^) tUj. The uniformizing coordinate is given by 
the potential by 

2 / 1 \ 3 , 



X = 63-(r , r „ 363^. (68) 



This is equivalent to (16 1 if one rescales into 6(j). The skeleton is a real half axis ImX — 0, ReAT < 63. 
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2. At the critical time t = 0, all roots coincide 63 = 62 = ei = 0. The curve further degenerates to a cusp 

= —iX^. The skeleton touches the double point. 

3. After the critical time i > 0, we push the evolution by splitting the double point e-i^ e\. The new branch points 
appear in the fluid, giving rise to shocks. By the previous argument, we only need to determine the constant 
term = 4e3|e2p. 

It follows that the roots and X scale as t^l"^ and Y scales as t^l'^. This elliptic curve has no deformation 
parameters, it is self-similar and in this sense unique. Scaling properties of the curve and the hydrodynamics 
are summarized as: 

X{u,t) = \t\^l-'Xr{\t\^l^u) , Yiu,t) = |t|3/4yi(|t|i/4„) ^ (69) 
Y{X,t)^\t\y^Y,(\t\-^/'x), ^{X,t) = \t\^^Ui{\t\-'/^X), (70) 

n{x, t) = \tf/*n, (itr^/'x) , vix, t) - \t\-'/*v, (itr^/'x) (71) 

where the functions subscripted by 1 are those evaluated at t = ±1. 

From the scaling properties alone and from the fact that Y^ is a polynomial of degree three it follows that 

2i 

n{X)^-i YdX^ {XY-2t(l3), (72) 

v{X) = G^^^^ , (73) 

where x| = ||f is a capacity (see below). We notice that after a critical time velocity diverges at branch points, 
where Y = 0. Eq. (72 1 allows to express the total mass deficit accumulated by a shock through the value of 
stream function at shock endpoints. Using ipi^es) = and Y(ei) = Y^e^) = 0, we have: 

j£'a|dX| = ^|Re0(ei)| = ^|V(ei)|. (74) 

The mass scales as t^^^. It grows faster than how liquid is getting drained (~ t). We will later show that the 
rate is a universal constant 

- a\dX\ « -(6.34513) t^/^ (75) 



B. Evolution of the elliptic curve 



The Krichever-Boutroux condition (61 1 uniquely determines and therefore evolution of the curve. In this 



case it requires that the integral over b-cycle be purely imaginary, 

ri(ei) = imaginary. (76) 



It follows from (72 1 that this is equivalent with vanishing pressure at branch points: 

p(ei,2) -Im0(ei,2) = 0. (77) 

Remarkably, but not accidentally, exactly the same problem appeared in a semiclassical analysis of triply- 
truncated solution of Painleve I equation [3H1 HQ] . 

To get the solution, let us parametrize the curve Y{X) by a uniformizing coordinate u as 

X = piu), (78) 
Y = ip'{u), (79) 

where p is the Weierstrass elliptic function whose half-periods Lo{t) and uj'{t) are (yet to be determined) complex 
functions of time. Since 63 is real, uj + lo' is real, and 62 = ei, so u> = Uj' . They scale with time as w ~ t~^^'^ . The 
branch points are given by (es, 62, ei) = (p(w -I- w'), p(u;), p(a;')). The rhombus-shaped fundamental domain is 
depicted in Figure [s] [4?. As t goes to zero, the rhombus becomes infinitely large, while preserving its aspect 
ratio. 
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FIG. 8: Fundamental domain (w-plane). The pressure level-lines (contour lines) and the stream lines (dashed) are drawn. 
Pre-image of shaded regions is a physical plane. Anti-Stokes lines (zero-level lines of Ref2) emanate from the branch points. 
Among them dashed (orange and green) lines are non-admissible anti-Stokes lines; the sign of Re Q is different on different sides 
of these lines. The thick (blue) lines are the admissible anti-Stokes lines - shocks. Refi is positive on both sides of these lines. 



The potential of the flow is obtained from the Darcy's law by 

HX) = -J^ dXdtY ^6i(^C{u)-^^^y}j , g2 = -12i. (80) 

Stream function, "0 = Hecj), and pressure, p — Im0, will be read from (t>{X). Here we have used the facts: 

^ dX , XdX , , 

C(«) = -i / T77^ > (81) 



Y{X) ' ' J Y{X) ' 



up to constants of integration. 



Eqs. (77 80 1 give the defining equation: 

3^3 _ + 
2 g2 uj + uj' 

whose solution is summarized as follows: 



J = u, (82) 



(61,62,63) ^ J ,1'^ (h+ih, i -i/i,-l)v^ w (0.276797-1- 1.797181,0.276797- 1.797181,-0. 553594)\/t, 
V 4ft^ — 3 ^ ^ 

92 = -l2t, 



53 = -l^y j' ^g + ^ t^/^ « -7.32176243U3/^ (83) 

where h « 3.246382253744278875676. This number comes from m = ^ + | , „ where m is the solution of 
the equation 

16m^ — 16m + 1 K{m) 



iim? — 9m + 1 Eim) 



(84) 



18 



that follows from (82 1. Here E and K are complete elliptic integrals. They also determine the invariant shape 

Imw K' {m) 



of the fundamental domain 



Rec 



K{m) 



0.81736372 



Eq. (84 1 follows from (82 1 with a help of formulas on p. 649 of 

K{m) 



92 



4(16m2 - lQm+l)K'^{m) 
3(w + 



- [6£'(m) + (4m - h)K{m)\ 



93 



8(2m - l)(32m2 - 32m - l)K'^{m) 



(85) 
(86) 



C. Shocks and anti-Stokes lines 



m 



II 



IV 



V 



FIG. 9: Geometry of anti-Stokes lines. The boundary condition Refi > on the remote part of the finger argX = tt, \X\ oo 
uniquely determines configuration of shocks r3,r2. 



Shock lines are level lines given by conditions (60 61 ) 



(87) 



Condition (61 1 



Im (XY) ^ 2iIm0(X) ^ 2tp{X), 

determines a total of seven anti-Stokes lines connected at the branch points. They are transcendental, computed 
numerically and denoted by Fi, Fy in Figure |9] Condition (88) is valid on both sides of a shock. Since Y changes 
sign through a shock, so does pressure. 

Among the seven anti-Stokes lines, only two F3, F2 are shocks, while Fi is a skeleton. They are selected by condition 
(6OI. Here is how it works. A remote part of the finger argX = tt, \X\ — > 00 is not affected by the critical transition. 
It is a skeleton. Therefore the branch of '--^ must be chosen such that Refi > in both sectors {IV, III}, 

where at large X, |argX — 7r| < 27r/5. This insures that a line density on the skeleton is positive. It follows that sign 
of Re H. is 



on {I, III, IV} 



on {II, V} 



(89) 



As a consequence, the signs are opposite on both sides of F4,F5,Fg and F7, and the same (plus) on both sides of 
ri,r2 and F3. The cuts selected this way satisfy the admissibility condition (60 1. 



Figure 10 illustrates the flow. There, level lines of Refi are drawn before, at, and after the critical time. 
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FIG. 10: Level lines of Re Q before, at, and after the critical time. The skeleton and shock lines are marked by (orange) lines 
(skeleton is hidden under the deep canyon), the branch points are marked by dots. For visually optimal presentation we have 
plotted arctan(Ref2). 



The Rankine-Hugoniot conditions (51 and (88 1 give the velocity of shocks. Noting that ImXY 
is a projection of a vector-coordinate oia point of a shock to the shock, we get 



(jXw , where X\\ 



(90) 



D. Capacity: discontinuous change of power 



The genus transition is signaled by an abrupt change of 53 . It follows from ( 83 ) that 

0. = Itp/^ I * (91) 

^3 I I \ -7.321762431, t > 0. ^ ' 

This discontinuity is related to a capacity C of the finger which is defined by the asymptote of the potential 4> at 
X ^00 (pTl): 



C 

^2x1/2 

For an elliptic curve, the capacity C is given by 



^{X) = -6iX'^' + i-^ + ... (92) 



3,„i/2/-8, t<0, 



C--g3- ^\t\ i 7.321762431, t > 0. ^^^^ 
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FIG. 11: Capacity passing the transition (continuous line). The dashed line corresponds to no discontinuity at the transition. 
The power N{t) — Nc = —Q^/(2nK) logC, features the same discontinuity. 



At a branching point, the capacity rate C/|t|^^^ goes through a discontinuous jump. Capacity passing through a 

Conformal radius of the fluid and power N(t) also exhibit the discontinuity (see 



critical point is drawn on Figure 
Sec. 
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We conjecture that the ratio of the time derivatives of the capacity before and after 1^2 branching, 

Caft cr branching 



7] :— lim 



cforc branching 



0.91522030388 



(94) 



is universal. It does not depend on the details of the flow. A discontinuity occurs every time when shocks branch 
This number can be used to express other objects of the curve. Let T] := f f at t > be: 



77 = sinh Gc, 0.820137. 



Then, for example, the branch points read: 

(ei,e2,e3) 



„ . , ©c + Tri _ . , 0c — Tri „ . , ©c 

2 smh , 2 smh , —2 smh — 

3 3 ' 3 



Vt. 



E. Detailed description of the shock 




FIG. 12: The close-up of shocks (left) and the line density profile of the shock a (right). In the upper panel, the shaded 
contours are equi-pressure lines and the dashed lines are the stream lines. The darker the shade is, the higher is the pressure in 
the fluid. The yellow dotted line is the zero-pressure line; it crosses the boundary of the fluid and a skeleton. With respect to 
this line, the main cut releases/absorbs fluid as a changes the sign. The arrows are moving directions of the shock and of the 
branch point. At 63 three directions: zero-pressure line, shock, and the velocity, are all at different angles; see the text. On the 
left of the zero-pressure line the finger (x, Y{x)) expands pushing the fluid away, on the right the finger retreats. In the right 
panel, we plot the line density a = |y(X)| for one of the shock (orange line in a vertical plane). The line density vanishes at 
branch points as a square root. The total charge on the shock grows like | (6.34513) i^^**. The blue line in the horizontal plane 
is the line density of the skeleton, and also a boundary of the fluid - a viscous finger. 



It is interesting to follow the zero-pressure line. It always emanates from the shock end points. Before the critical 
time, it is also the boundary of the fluid - a flnger. It is given by the graph {x,Y{x)), where x — KeX. After the 
critical time, a zero pressure line has two branches. One remains the boundary of the fluid, while another branch is 
in the fluid and crosses the boundary at a point where dcj) = 0, as on the upper panel of Figure [T2} This occurs at 
X — — ^Tjq^^^ ~ ~if| ~ —rjt^^'^. The zero-pressure line emanates from the branch point ei. If ei is choosen as origin, 
the zero-pressure line emanates with the angle ^zcro-prcssuro ~ 0.2357r, as we will evaluate below. 

The significance of the point d(/) = and the zero-pressure line can be observed in Figure [12] One can see that the 
stream lines (dotted lines) of fluid that start from the real cut (a skeleton) , flow in the opposite directions with respect 
to the zero-pressure line: the one on the right goes toward the real cut, the one on the left goes away. Exactly at 
d0 = 0, the flow lines are intersecting perpendicularly. This means that the boundary of the fluid (a graph {x, Y{x))) 
on the left of this point retreats toward the drain, while the finger expands, as it was before the critical time. From 
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the right to this point, the boundary moves in the opposite direction, the finger retreats smoothing the tip. The point 
itself is a stationary point of fluid where the velocity of flow and therefore ViidiscT/" vanishes. It follows from (56 1 
that a also vanishes at this point. 

As has been pointed out, the evolution of the curve Y(X') is self- similar: the whole picture in Figure |7] simply 
expands with respect to the origin by the factor \/t for t > (or for i < 0). Therefore, the branching shocks are 
moving to the darker region (to the left) as time goes on (see the big arrow for the motion of the shock). 

In Figure |7] anti-Stokes lines are located by the condition Re = 0. The bold lines are admissible anti-Stokes 
lines selected by the condition Re > on both sides of the fluid. We have argued that this condition leads to the 
curl-free flow of the entire fluid which, after the critical time, consists of decompressed fluid in shocks and a fluid with 
a constant density elsewhere. 

Let us discuss the behavior of physical quantities around branch points, say ei. The line density is 

a^\Y\^2\j\V~W\X -ei\ , (95) 
where I7I « 2.66757 comes from the estimate below: 

7= -j:V{e2 - ei)(ei - 63) « 2.60534- 10.57281. 

From this we can estimate the angle of the shock line ^shock at its end point ei, by 

^fshock + arg 7 = I e.hock ~ 0.04592497r + ^ . 



Here and below all angles are with respect to the real axis 



The angle between two shocks at their ends is twice (0.09184987r+ ^). It is by 16.532964° larger that the ^ angle 



between shocks at their origin. 

We also estimate the density near 63 on a skeleton x < 63. It is (7(2:) = |y| 2 •\/— (63 — ei)(e3 — 62) ^/e^ — x. It 
changes in a discontinuous manner similar to the capacity, ( 93 1 

lim ^^fter branching ^ 1^3 - ^ ^^^^^^^ ^^^^ 
x^es a^Q^Qj-Q branching 3 v i 



The graph of density is depicted on the lower panel in Figure [12] 



The total mass deficit carried by a shock follows from ( 74 1 and ( 80 ) 



/ .|dXH-6 



^(6.34513) i^/^ 




Simple calculations give the potential (pressure and stream function) at the branch point 

(j)iX) - ei-/ty/X - ei « (1.7506 + 14.5237) ty/X - ei . 
This gives the angle of the zero-pressure line at ei 

^6'zcro-prossuro + arg(ei7) = ^ 6'zoro-prcssure = 0.23506l7r , (97) 

such that the angle between the two zero-pressure lines that emanate, respectively, from ei and 62, is twice the above: 
0.4701227r, slightly less than n/2. 

Finally, an angle of the velocity of the branch point ei (determined by the scaling relation (69 ) efe = ^e^;, fc = 1, 2, 3) 
is Oyeiocity = 0.4513577r. 

Summing up, velocity of an end point, zero-pressure line and the shock itself at the end point ei (and symmetrically 
at 62) all have different directions. The relative angles between them are 0^ — Ozao pressure = 0.2162967r, 9y — ^?shock = 
0.07209877r. The velocity of the end point is depicted by an arrow in Figure [T2| 



VII. DISCUSSION 



For most initial configurations, zero-surface Hele-Shaw flows evolve into a cusp-singularity in finite time. As such 
the problem is ill-defined. A singularity signals that the zero-surface Hele-Shaw flow is a singular limit of a realistic 
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and a well-defined problem, where velocity and pressure rapidly change at a small scale controlled by some additional 
small parameters. In this limit, when these parameters are set to zero, areas of rapidly changing gradients shrink to 
lines - shocks. There, velocity and pressure are not differentiable functions and may even diverge. 

In this paper we found a weak solution of the Hele-Shaw problem, where a few simple physical assumptions almost 
uniquely determine the evolution of the flow through a singularity, to a growing and branching shock graph. 

Hele-Shaw flow occurs in a variety of hydrodynamics settings [m [161 E]) aggregation models [18] and also in 
quantum electronic liquids [43 , random matrices, theory of orthogonal polynomials [33J and complex analysis [441 145] . 
In all cases except hydrodynamics a small parameter is available: in aggregation problems it is a size of a aggregating 
particle, in quantum problems a small parameter is h, in problems with random matrices and orthogonal polynomial 
a small parameter is 1/N, where TV is a size of the matrix or an order of a polynomial. In hydrodynamics, ithe 
varies and depends on a setting. The latter gives a particular (setting-dependent) interpretation of the viscous shock 
solutions and requires a case by case study. 

In the problem with two fluids (an inviscid fluid pushing a viscous fluid), a shock is a narrow channel filled by a 
compressed inviscid fluid which forms a line of vortices. In the setting with one viscous fluid (such as a thin layer 
on a wet surface), sucked away at a distant point, shocks are narrow channels where a fluid layer cracks such as the 
wet surface (substrate) supplies a turbulent and compressed fluid to the layer. In both cases, one must relax the 
condition of incompressibility and zero-vorticity flow. Fluids are compressible and not curl-free at the shock scale. 
The major physical principle which determines the weak solution is the requirement that the entire fluid (or two 
fluids) be curl-free and incompressible at a large scale. This means that the smooth part of the flow adjusts to keep 
the total vorticity zero and a constant draining rate. 

In either case, a well-deflned set of hydrodynamics equations which may lead to a shock solution is lacking. 

Contrary to a differential Darcy law, the weak solution is not formulated as a Cauchy problem. It relays on an 
assumption that the boundary fluid is an algebraic curve. In this case, the weak solution uniquely determines the 
evolution. Although it seems a rather limited set of initial conditions, the algebraic curve occurs locally at the cusp 
singularity. In this manner the weak solution describes the evolution through a singularity. 

The condition that the entire flow (a smooth part of the flow and shocks combined) is curl-free brings the problem to 
an evolution of Krichever-Boutroux curve. These are very special curves which previously appeared in a semiclassical 
analysis of certain solutions of Painleve I equation [3HJ [101 [lH HI] • 

In the paper we studied the most generic (2,3) cusp singularity where the Krichever-Boutroux curve is elliptic. 
In this case the flow produces two shocks. This solution is universal (parameter-independent) and self-similar. We 
emphasize interesting universal numbers which describe jumps of physical variables (like the power N(t) and capacity 
C{t) (93) when the flow goes through a singularity. 

Our solution represents a local branching event of a further developed tree. In Figure [2] we depicted a numerical 
solution with two generation of branching. An unknown, interesting global structure of the shock's branching tree 
is far out of the scope of this paper. It will be very interesting to see whether a developed shock's tree exhibit a 
universal scale invariant limit after a large number of branchings. 
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